Cite this Article
Zhao Lei, Zhang Cun-bo, Liu Jian-xin, Luo Ji-sheng. Improved algorithm for solving nonlinear parabolized stability equations. Chinese Physics B, 2016, 25(8): 084701  
Permissions
Improved algorithm for solving nonlinear parabolized stability equations
Zhao Lei, Zhang Cun-bo, Liu Jian-xin†, , Luo Ji-sheng
Department of Mechanics, Tianjin University, Tianjin 300072, China

 

† Corresponding author. E-mail: shookware@tju.edu.cn

Project supported by the National Natural Science Foundation of China (Grant Nos. 11332007 and 11402167).

Abstract
Abstract

Due to its high computational efficiency and ability to consider nonparallel and nonlinear effects, nonlinear parabolized stability equations (NPSE) approach has been widely used to study the stability and transition mechanisms. However, it often diverges in hypersonic boundary layers when the amplitude of disturbance reaches a certain level. In this study, an improved algorithm for solving NPSE is developed. In this algorithm, the mean flow distortion is included into the linear operator instead of into the nonlinear forcing terms in NPSE. An under-relaxation factor for computing the nonlinear terms is introduced during the iteration process to guarantee the robustness of the algorithm. Two case studies, the nonlinear development of stationary crossflow vortices and the fundamental resonance of the second mode disturbance in hypersonic boundary layers, are presented to validate the proposed algorithm for NPSE. Results from direct numerical simulation (DNS) are regarded as the baseline for comparison. Good agreement can be found between the proposed algorithm and DNS, which indicates the great potential of the proposed method on studying the crossflow and streamwise instability in hypersonic boundary layers.

PACS: 47.20.Ib;47.40.Ki;47.20.Lz
Keyword:improved NPSE approach;stationary crossflow vortices;fundamental resonance;hypersonic;boundary layers
1. Introduction

The mechanisms of transition from laminar to turbulent flow are of strong academic interests. Understanding the transition process and predicting its location are necessary for accurate computing of aerodynamic drag and heat protection configurations. Under flight conditions, natural transition is dominant because of the low free stream turbulence level. Therefore, stability theory plays a significant role in understanding the transition. There are usually three methods to describe the amplification of unstable disturbances in boundary layers. They are linear stability theory (LST), direct numerical simulation (DNS), and (nonlinear) parabolized stability equations (PSE/NPSE). LST is a local method based on the quasi-parallel hypothesis. It can predict the linear growth of an infinitesimal disturbance well.[1,2] However, it ignores the nonparallel and nonlinear effects. DNS can obtain complete information of disturbances through solving the Navier–Stokes equations directly but with high computational demands. NPSE approach developed by Bertolotti and Herbert[3,4] takes into account the nonparallel and nonlinear effects, the solution accuracy of NPSE is comparable with DNS while the computational cost is much less.

NPSE is solved by a marching method. Due to its high computational efficiency and accuracy, The NPSE approach has been widely used to study the stability and transition. This approach was initially applied to an incompressible Blasius boundary layer by Bertolotti.[3,4] Then its application was extended to three-dimensional flow and compressible flow.[5] For the crossflow instability, Wang et al.[6] applied this approach to calculate the development of crossflow waves and predicted the nonlinear saturation characteristic. Haynes[7] made detailed comparisons between NPSE results and experimental measurements.[8] The remarkably good agreements validated this approach for three-dimensional flow subjected to crossflow instability. Malik[9] adopted this approach to simulate the nonlinear development of stationary crossflow vortices over a swept airfoil, and investigated the secondary instability of saturated crossflow vortices. Xu[10,11] and Li[12,13] also obtained the baseflow for the secondary instability analysis by this approach.

For investigating the transition dominated by streamwise instabilities, NPSE approach is also an efficient method. By this approach, Zhang and Zhou[14] investigated the subharmonic resonance in a subsonic boundary layer and acquired an essential agreement with DNS results.[15] Mayer[16] simulated the oblique breakdown initialized by first-mode disturbances in a Ma = 3 boundary layer by NPSE approach. The transition onset location agrees well with DNS results.[17] Zhang and Zhou[18] also investigated the subharmonic secondary instability induced by first-mode and second-mode disturbances in a Ma = 4.5 boundary layer by this approach.

Despite the computational advantage of the NPSE approach, numerical instabilities may lead to convergence failure when the disturbance amplitude increases to a threshold,[19] especially in hypersonic flow.[20] Chang[21] and Li[12] calculated the development of stationary crossflow vortices in a Ma = 2.4 supersonic flow over a swept airfoil. For the most unstable mode disturbances, all calculations diverge when the amplitudes of the disturbances reach the levels of about 0.17. Kuehl[22] also investigated the development of stationary crossflow vortices in a hypersonic boundary on a cone with 6 degree yaw. The NPSE calculation also diverges in downstream when the amplitude reaches about 0.13. The work of Kuehl concludes that the NPSE approach tends to diverge when the amplitudes of stationary crossflow disturbances increase to a certain level in hypersonic flow. However, the saturated stationary crossflow vortices are required to be considered in order to study the breakdown induced by the secondary instability as in the case of low-speed or/and subsonic flow.

For this aim, therefore, an improved algorithm for solving NPSE is developed in this research. The nonlinear development of stationary crossflow vortices and the fundamental resonance of second-mode disturbances are investigated in hypersonic boundary layers. The results obtained by the improved NPSE approach were compared with those by the traditional NPSE approach and DNS for validation.

2. Methodology
2.1. Traditional NPSE

In general, the instantaneous flow quantities q = [ρ,u,v,w,T] can be decomposed into steady laminar base flow and perturbation quantities in the form q = q0 + q′. Here, ρ is the density, T is temperature, u,v,w represent the velocity components in streamwise, wall-normal, and spanwise directions respectively; q0 = [ρ0,u0,v0,w0,T0] are the base flow quantities and q′ = [ρ′,u′,v′,w′,T′] the corresponding perturbation quantities. The disturbance equations can be derived from Navier–Stokes equations. After a nondimensional procedure they are presented in a compact form as follows:

where x,y,z represent streamwise, wall-normal, and spanwise directions respectively and t denotes the time. Coefficient matrixes Γ, A, B, C, D, Vyy, Vxx, Vzz, Vxy, Vxz, Vyz are functions of base flow quantities. Nonlinear terms are denoted by the vector F which is a function of perturbation and base flow quantities.

Each perturbation quantity is assumed to be periodic in temporal and spanwise directions; the disturbances and nonlinear terms are then expressed by the truncated Fourier series

where M1 and M2 represent one-half of modes used in the analysis; β is the nondimensional spanwise wave number, and ω the nondimensional frequency of a disturbance. is the shape function vector of the Fourier mode (m,n) and αmn the streamwise complex wave number; they both vary slowly in x direction. Substituting Eqs. (2a) and (2b) into Eq. (1), the stability equations are derived as

The matrixes Amn, Bmn, and are given by

By analyzing the magnitude of each term in Eq. (3), the partial derivative of viscous terms at x direction, i.e., and can be ignored. In this manner, the nonlinear parabolized stability equations are obtained.

This equation has two unknown variables, and αmn, therefore an additional condition is required to remove the arbitrariness. This supplementary condition is usually not unique. In three-dimension compressible boundary layer the following form is applied,

In the solution procedure, equation (6) is utilized in the form of iterative formula below,

where

the superscript * denotes the complex conjugate.

Equation (5) can be rapidly solved by a marching method. Discretizing the partial derivative will eventually lead to a system of algebraic equations as

in which the superscript i + 1 denotes the current streamwise location where the perturbation quantities are to be solved and the i denotes the last streamwise location where all perturbation quantities are known. The matrixes Lmn and Rmn consist of the matrixes appearing in Eq. (5) and the finite difference operators. They are both unary functions of q0. The vector is a binary function of q0 and the unknown variable q′. Therefore, an iterative procedure for is needed to solve the

The marching approach for solving Eq. (8) can be outlined as follows:

2.2. Improved NPSE approach

In hypersonic boundary layers, the traditional NPSE approach easily diverges when the nonlinear interactions are significant. Chang[19] pointed out that these significant nonlinear interactions lead the nonlinear forcing terms to carry more weight than linear operator in the NPSE. Dong[23] and Cao[24] found that the mean flow distortion has a significant contribution to the instabilities of higher harmonic modes. In their work, the mean flow distortion was added to the base flow and the stability of the distorted base flow was study. The results from their work indicate the superiority of adding the flow distortion in a linear operator rather than the nonlinear forcing terms. What is more, moving the mean flow distortion into the linear operator benefits to reduce the weight of the nonlinear forcing terms. The improved algorithm for NPSE proposed by this paper is based on this idea.

Different from the traditional approach, the improved algorithm considers the mean flow distortion in the linear operator. The iterative form governing equations of the improved algorithm can be represented as

Here, is the mean flow distortion induced by the nonlinear interactions. In this algorithm, the mode (0,0) is solved by the same method with the traditional NPSE approach. However, the mean flow distortion is moved into the linear operator instead of the nonlinear terms when solving the higher harmonic modes. In other words, the base flow and the mean flow distortion are both contained in the linear operator and have linear contributions to the other modes.

Additionally, an under-relaxation factor is introduced in the formula (4c) for robustness when calculating iteratively. The new form of the formula (4c) containing an under-relaxation factor can be written as

where Fmn(old) and Fmn are the values obtained by the formula (2b) in the previous and current iterations; δ is the under-relaxation factor set as 0.2.

3. Results and discussion

The improved NPSE approach is applied to two challenging research cases in this study. One is the nonlinear development of stationary crossflow vortices in a hypersonic three-dimensional boundary layer; the other is the fundamental resonance of second-mode disturbances in a hypersonic flat-plate boundary layer. In order to verify the effectiveness and reliability of this improved approach, the results by traditional NPSE approach and by DNS are also provided as comparisons.

3.1. Nonlinear development of stationary crossflow vortices

Accurately simulating the nonlinear developments and amplitude saturation characteristics of stationary crossflow vortices is a necessary work for studying the crossflow secondary instability. In hypersonic three-dimensional boundary layers, however, the traditional NPSE algorithm tends to diverge when the amplitude of disturbance amplifies to a threshold.[21,22]

3.1.1. Steady base flow and linear stability analysis

Figure 1 sketches the physical model, a swept blunt plate with infinite spanwise length is yawed in hypersonic free steam. The temperature of the free stream T∞ = 226.5 K, the Mach number Ma = 10, and the Reynolds number Re = ρ∞Q∞R0/μ∞ = 131800, where R0 = 35 mm is the nose radius at the leading edge of the blunt plate. The sweep angle is SA = 45°. The flow quantities are scaled by the free stream values. ρ = ρ*/ρ∞, u = u*/Q∞, v = v*/Q∞, w = w*/Q∞, T = T*/T∞, The coordinates are scaled by the nose radius R0, x = x*/R0, y = y*/R0, z = z*/R0. The hypersonic steady flow over the swept blunt plate is obtained by an NS equations solver. The hypothesis of calorically perfect gas is applied. An adiabatic and no slip condition is adopted on the wall.

Fig. 1. Schematic illustration of the physical model.

Figure 2 displays the development of streamline and crossflow velocity components at different streamwise locations. It shows that the boundary layer thickness increases in the streamwise direction while the crossflow velocity decreases slowly. The crossflow instability is analyzed by LST, and two unstable disturbances with different spanwise wavelengths are investigated. Their spanwise wavelengths are 109.96 mm and 54.98 mm, respectively.

Fig. 2. Velocity distribution at x*/R0 = 6, 12, 18, 24, 30, 36, 42 (The arrow marks the increase of x* location): (a) Streamline velocity profiles; (b) Crossflow velocity profiles.

Figure 3 shows their growth rates and N factors. As shown in the figure, the disturbance with λz = 109.96 mm has a large growth rate in the whole computational zone, and its N factor increases to 10 downstream. The disturbance with λz = 54.98 mm amplifies by a larger growth rate upstream, and then it damps after the location of x = 57 where its N factor reaches a peak value of 7. These two unstable disturbances have different growth characteristics, and their large N factors mean that the disturbances can amplify to a saturated amplitude. The traditional NPSE approach will diverge in these cases. These are the reasons for setting the physical parameters as described above.

Fig. 3. The growth rates and N factors of two waves obtained by LST: (a) growth rates and (b) N factors.
3.1.2. Nonlinear development of stationary crossflow vortices

The inlet conditions for NPSE and DNS are provided by the LST results. For the unstable disturbance with the spanwise wavelength 54.98 mm, the N factor (see Fig. (3b)) suggests the disturbance will amplify to a peak value at about x = 57. The nonlinear evolution of disturbances with two different initial amplitudes are computed. The initial amplitudes of disturbances (based on the streamwise velocity disturbance) are specified as A0 = 8×10−5 and A0 = 1×10−3 respectively.

Figure 4 displays the amplitude evolutions of different modes. For the case of A0 = 8×10−5, as can be seen in Fig. 4(a), the amplitude of the fundamental mode (0,1) increases to a peak value of about 0.07 at about x = 48, and then it decreases quickly. The damping location of the disturbance is brought forward due to the nonlinear interactions. Because the amplitude is not large, both of the traditional NPSE and improved NPSE approaches accurately simulate the development of amplitudes and the agreement with the DNS result is perfect.

Fig. 4. The nonlinear evolution of the amplitudes of streamwise velocity disturbance (λz = 54.98 mm) for (a) A0 = 8×10−5 and (b) A0 = 1×10−3. The vertical line denotes the divergence location of the traditional NPSE approach.

For the case of A0 = 1×10−3, as shown in Fig. 4(b), the amplitude of the fundamental mode (0,1) reaches a peak value of about 0.21 at x = 32, and then decreases quickly. The mean flow distortion reaches a peak amplitude of about 0.12. In this case, the disturbance amplifies to a larger amplitude compared with the small initial value case and the nonlinear interactions are also more significant. It can be seen that the result of the improved NPSE approach is in excellent agreement with that of DNS. However the traditional NPSE approach diverges at the location of x = 26 where the amplitude of the mode (0,1) just reaches to 0.13.

For the unstable disturbance with the spanwise wavelength 109.96 mm, the LST result (see Fig. 3(b)) suggests the disturbance amplifies in the whole region. As in the case above, two cases with different initial amplitudes are also investigated. The initial amplitudes are A0 = 3×10−5 and A0 = 1×10−3 respectively. The amplitudes evolution of the streamwise velocity disturbances are presented in Fig. 5. For the case of A0 = 3×10−5, the disturbance amplifies rapidly in a long streamwise range, and the amplification tends to slow finally due to the strong nonlinear interactions. The amplitude of the mode (0,1) reaches to around 0.25 finally. As shown in Fig. 5(a), the result obtained by the improved NPSE approach agrees excellently with that of DNS. However the traditional NPSE calculation diverges when the amplitude of the mode (0,1) increases to 0.11 at x = 69.

Fig. 5. The nonlinear evolutions of the amplitudes of streamwise velocity disturbance (λz = 109.96 mm) (a) A0 = 3×10−5 and (b) A0 = 1×10−3. The vertical line denotes the divergence location of the traditional NPSE approach.

When the initial amplitude increases to A0 = 1×10−3, as shown in Fig. 5(b), the disturbance amplifies quickly, and then the amplitude of the mode (0,1) sustains equilibrium at a level of about 0.23. The improved NPSE approach still obtains an excellent result comparable to DNS except for a slight deviation in the tail of the computational field. However, the traditional NPSE approach diverges at the location x = 36 where the amplitude of the mode (0,1) just reaches 0.13.

As shown in Figs. 4 and 5, the traditional NPSE approach diverges easily when the amplitude of the mode (0,1) exceeds 0.1. It is also reported in Kuehl’s work[22] when calculating the nonlinear development of stationary crossflow vortices in hypersonic boundary layers. However, it is shown that the improved NPSE approach still excellently simulates the development of disturbances even though the amplitude of the mode (0,1) increases to about 0.25. The amplitude saturation characteristic of stationary crossflow vortices can also be obtained by this approach. To obtain the results presented in Figs. 4 and 5, the improved NPSE approach took no more than 15 minutes CPU time however DNS required more than 400 hours.

To further examine the improved NPSE approach, the eigenfunctions development of streamwise velocity disturbance are studied. In Fig. 6, DNS results are used as baseline for comparison. As shown in which, the improved NPSE approach accurately describes the evolution of the fundamental mode (0,1), the mean flow distortion mode (0,0), and the harmonic mode (0,2) along the streamwise direction. The results of the improved NPSE agree excellently with those of DNS.

Fig. 6. Development of the streamwise velocity disturbance in x direction (λz = 109.96 mm, A0 = 3×10−5). The arrow marks the increase of streamwise location x.
3.2. Fundamental resonance of second mode disturbances

When the traditional NPSE approach is applied to calculate the fundamental resonance of second-mode disturbances, it usually diverges before the 3D disturbance amplifies to a large amplitude. Therefore, the computational work on this subject is mainly done by DNS rather than by NPSE approach. In this section, the improved NPSE approach is applied to study the fundamental resonance of second mode disturbances in a hypersonic boundary layer.

Yu[25] investigated the nonlinear interaction mechanisms of disturbances in supersonic flat-plate boundary layers by DNS. In Ref. [25], the case1 in the table 1 illustrated the nonlinear interaction mechanisms of Klebanoff-type second mode disturbances. The same case is investigated by the improved NPSE approach in this study. An unstable 2D disturbance of the second mode is introduced at the inlet, accompanied by a pair of stable 3D disturbances with the same frequency. Figure 7 shows the evolutions of Klebanoff-type second mode disturbances in a Ma = 6 flat-plate boundary layer. As illustrated by the DNS results, the 3D disturbance is damped at the beginning as predicted by the linear theory. When the 2D disturbance increases to a certain level, the 3D disturbance amplifies rapidly from x = 550 due to the effect of fundamental resonance. The amplitude of the 2D disturbance at this location is about 0.09. As shown in the figure, the traditional NPSE approach diverges at x = 558 where the 3D disturbance just starts to amplify. However the improved NPSE approach describes the process of fundamental resonance excellently. Despite its success in studying the dynamics of the fundamental resonance, we should point out that the improved algorithm diverges downstream after x = 590.

Fig. 7. The evolution of Klebanoff-type second mode disturbances in Ma = 6 flat-plate boundary layer. The vertical line denotes the divergence location of the traditional NPSE approach.
4. Conclusions

An improved algorithm for solving NPSE is developed in this paper to address the divergence problem of the traditional NPSE solving technique for hypersonic flow. In this approach, the mean flow distortion generated by nonlinear interaction is included into the linear operator. An under-relaxation factor is introduced in the iterative computation of the nonlinear terms. This improved NPSE approach is applied to study the nonlinear evolution of disturbances induced by the crossflow instability and streamwise instability in hypersonic boundary layers. Good agreement with DNS is achieved in the challenging cases. The improved algorithm is able to accurately capture the large-amplitude saturation characteristic of stationary crossflow vortices and the process of fundamental resonance in hypersonic flow. It is shown that the proposed algorithm has great potential in studying the crossflow instability and streamwise instability in hypersonic boundary layers.

Reference
1Mack L M1984“Spec. Course Stab. Transit. Laminar Flows”AGARD Rep.709pp. 3.1–81
2Reed H LSaric W SArnal D 1996 Ann. Rev. Fluid Mech. 28 389
3Bertolotti F P1991“Linear and nonlinear stability of boundary layers with streamwise varying properties” Ph. D. DissertationColumbusOhio State University
4Bertolotti F PHerbert TSpalart P R 1992 J. Fluid Mech. 242 441
5Herbert T 1997 Ann. Rev. Fluid Mech. 29 245
6Wang M JHerbert TStuckert G199425th AIAA Fluid Dynamics Conference1994Colorado Springs, USA10.2514/6.1994-2374
7Haynes T SReed H L 2000 J. Fluid Mech. 405 325
8Reibert M SSaric W SCarrillo R BChapman K L199634th AIAA Aerospace Sciences Meeting and ExhibitJanuary 15–18, 1996Reno, NV, USA18410.2514/6.1996-184
9Malik M RLi FChoudhari M MChang C L 1999 J. Fluid Mech. 399 85
10Xu G LXiao Z XFu S 2011 Sci. China-Phys. Mech. Astron. 54 724
11Xu G LXiao Z XFu S 2011 Sci. China-Phys. Mech. Astron. 54 2040
12Li FChoudhari M M 2011 Theor. Comput. Fluid Dyn. 25 65
13Li FChoudhari M MDuan LChang C L 1994 Phys. Fluids 26 064104
14Zhang Y MZhou H 2008 Appl. Math. Mech.-Engl. Ed. 29 833
15Li N2007“Investigation on transition of laminar-turbulent in a boundary layer on a flat plate based on the spatial model” Ph. D. DissertationTianjinTianjin University(in Chinese)
16Mayer C S JFasel H FChoudhari M MChang C L 2014 AIAA 52 882
17Mayer C S JDominic AFasel H F 2011 J. Fluid Mech. 674 5
18Zhang Y MZhou H 2008 Appl. Math. Mech.-Engl. Ed. 29 1
19Chang C L2004NASA Langley Research Center, NASA TM 213233
20Ren JFu S 2015 J. Fluid Mech. 781 388
21Choudhari M MChang C LStreett CBalakumar P2003AIAA 41st Aerospace Sciences Meeting and ExhibitJanuary 6–9, 2003Reno, Nevada, USA97310.2514/6.2003-973
22Kuehl JPerez EReed H L201250th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace ExpositionJanuary 9–12, 2012Nashville, Tennessee, USA92110.2514/6.2012-921
23Dong MLuo J S 2007 Appl. Math. Mech.-Engl. Ed. 28 1019
24Cao WHuang Z F 2006 Appl. Math. Mech.-Engl. Ed. 27 425
25Yu MLuo J S 2014 Sci. China-Phys. Mech. Astron. 57 2141